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Abstract 

The long-term distributions of trajectories of a flow are described by invariant densities, 
i.e. fixed points of an associated transfer operator. In addition, global slowly mixing struc- 
tures, such as almost-invariant sets, which partition phase space into regions that are almost 
dynamically disconnected, can also be identified by certain eigenfunctions of this operator. 
Indeed, these structures are often hard to obtain by brute-force trajectory-based analyses. In 
a wide variety of applications, transfer operators have proven to be very efficient tools for an 
analysis of the global behavior of a dynamical system. 

The computationally most expensive step in the construction of an approximate transfer 
operator is the numerical integration of many short term trajectories. In this paper, we 
propose to directly work with the tnfinitesimal generator instead of the operator, completely 
avoiding trajectory integration. We propose two different discretization schemes; a cell based 
discretization and a spectral collocation approach. Convergence can be shown in certain 
circumstances. We demonstrate numerically that our approach is much more efficient than 
the operator approach, sometimes by several orders of magnitude. 

1 Introduction 

Analysis of the long-term behavior of flows can be broadly classified into geometric methods 
and statistical methods. Geometrical methods include the determination of fixed points, periodic 
orbits, and invariant manifolds. Invariant manifolds of fixed points or periodic orbits act as barriers 
to transport as trajectories may not cross the manifolds transversally. Statistical methods include 
determining the distribution of points in very long trajectories of a very large set of initial points 
(i.e. a physical invariant measure |29) . often possessing an invariant density) and the identification 
of meta-stable or almost-invariant sets [51 [TSl [T3] . Almost-invariant sets partition the phase space 
into almost dynamically disconnected regions and are important for revealing global dynamical 
structures that are often invisible to an analysis of trajectories. These metastable dynamics also 
go under the names of persistent patterns, or strange eigenmodes, both of which are realisable as 
eigenfunctions of a transfer operator. Frequently, the boundaries of maximally almost-invariant 
sets (those sets which are locally closest to invariant sets) coincide with certain invariant manifolds 

m- 

Our focus in the present paper is on statistical methods, although we demonstrate via case 
studies the relationships to geometrical methods. The most commonly used tool for statistical 
methods is the transfer operator (or Perron- Frobenius operator). Fixed points of the transfer 
operator correspond to invariant densities, while eigenfunctions corresponding to real positive 
eigenvalues strictly less than one provide information on almost-invariant sets. In practice, one 
typically constructs a finite-rank numerical approximation of a transfer operator and computes 
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large spectral values and eigenfunctions for this finite-rank operator. The construction of the 
finite-rank approximation requires the integration of many relatively short trajectories with initial 
points sampled over the domain of the flow. It is this use of short trajectories that gives the 
transfer operator approach additional stability and accuracy when compared with computations 
based upon very long trajectories. Long trajectories continually accumulate small errors from 
imperfect numerical integration and finite computer representation of numbers; these small errors 
quickly grow in chaotic flows. While the transfer operator approach is very stable, it still requires 
the computation of many small trajectories which can be very time consuming in some systems. 
The approach we describe in the present work obviates the need for any trajectory integration at 
all and works directly with the vector field. 

Our approach exploits the fact that the evolution of probability densities u = u(t) can be de- 
scribed by generalized solutions of the abstract Cauchy problem dtu — Au. The Perron-Frobenius 
operator is the evolution operator of this equation, and has the same eigenfunctions as the opera- 
tor A. The operator A is an unbounded hyperbolic (if the underlying dynamics is deterministic) 
or elliptic (if the deterministic dynamics is perturbed by white noise) partial differential operator. 
Standard techniques allow us to approximate the eigenmodes we are interested in: finite difference, 
finite volume, finite element and spectral methods yield such discretizations; see [251 113 EH HI [5] 
and the references therein. 

An outline of the paper is as follows. In Section 2 we provide background on the infinitesi- 
mal operator arising from smooth vector fields and describe conditions under which the operator 
generates a semigroup of transfer operators. In order to obtain formal results, we will require 
the addition of a small amount of diffusion to the deterministic flow. In the small diffusion set- 
ting we discuss existence of invariant densities and spectral results for the associated infinitesimal 
operator. Section 3 describes a spectral Galerkin method for the approximation of the infinitesi- 
mal operator. We apply results from the numerical analysis of advection-diffusion PDEs to show 
that our Galerkin method approximates the true eigenfunctions of the infinitesimal operator as 
our Galerkin basis becomes increasingly refined. We can also show that the convergence rate is 
spectral; that is, faster than any polynomial. Section 4 describes an Ulam-based Galerkin method 
for approximating the infinitesimal operator. This Ulam-based approach is new and shares some 
similarities with finite-difference schemes. While we cannot show convergence of this approxima- 
tion scheme, the numerical results obtained are extremely fast and accurate. Sections 5-8 detail 
the practical application of our two approximation methods to flows in one-, two-, and three- 
dimensional domains. We demonstrate the spectral accuracy of our spectral Galerkin method and 
compare with the accuracy of the Ulam-based Galerkin method. We also compare the accuracy vs. 
computational effort of our two new approaches with standard transfer operator approaches. The 
natural relationships between the outputs of our statistical methods and geometric objects such 
as the vector field and invariant manifolds are also elucidated in each case study. We conclude in 
Section 9. 

2 Dynamics, densities and semigroups 

Let the domain M C M'' of our flow be a smooth compact manifold and m the (normalized) 
Lebesgue measure on M . Denote hy F : M ^ the vector field generating the flow and by 
<I>* : M M, t £ R, the flow, i.e. <I>*(x) represents the location of a trajectory beginning at 
X G M after flowing for t G M time units. One has that d^*{x)/dt = F{x). Note that - provided 
that the components Fi,i ^ 1, . . . ,d have continuous derivatives ~ the function 4>* : M — >■ M is a 
diffeomorphism for every t g M. 

Invariant sets are structures of dynamical interest. A set A G M is called invariant, iff 
<I>^*^ = A for all t. Also, one asks how the flow changes probability measures. Sample x according 
to a probability measure /i; the distribution of <I>*a; is then given hy ^ o <I>~*. Special attention 
is to be drawn to invariant measures, which do not change under the dynamics {fj, = ^ o $~*). 
Invariant measures fi are called ergodic if invariant sets have either zero or full measure, i.e. if 
A C M satisfies = A then ii{A) £ {0, 1}. An even more restricted class of ergodic measures 
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are the physically relevant (or natural) ones, satisfying 



1 



L 



T 



"0 d/i = lim 



T 



dt 



for all ip : M ^ M. continuous and x E U C M with m{U) > 0. One has that absolutely 
continuous!^ ergodic measures are natural invariant measures. The density of an invariant measure 
is called the invariant density. 

2.1 Transfer operators 

When looking for invariant densities, one can rephrase the action of the flow on measures as an 
action on densities. If / denotes the density of /i, then / is evolved by the flow as 



where \B\ denotes | det B\ for a matrix B. The linear operator 7^* : L^{M) O is known as a transfer 
operator or the Perron- Frohenius operator associated with the flow $. Note that invariant densities 
are fixed points of 7-"* . For each t € IR we have 

(i) 7" is linear, 

(ii) 7^*/ > if / > and, 

(iii) \\V*'f\\ = 11/11 for ah / > 0, where |1 • |1 is the norm. 

i.e. is a Markov operator. Moreover, these properties also imply WV^W < 1, hence "P* is a 
contraction and its eigenvalues lie inside the unit disc. In fact, since the flow $* is a bijection, we 
have = ||/|| for all / G L^. To see this, note that 



for i > and A C M measurable. Now write f = — , with f^{x) = max{±/(a;), 0}, then 
11/11 — /(/^ — /^) = Jg+ /^ + Jg- /^, where S"^ is the support of /=*=. Now use properties (ii), 
(iii) and ([T]) with A = S^. A consequence of this is that the eigenvalues of V* lie on the unit 
circle. 

2.2 Operator semigroups, generators 

The transfer operator also inherits some (semi)group properties of the flow 

Definition 2.1. Let {X, || • ||) be a Banach space. A one parameter family {T*}f>Q of bounded 
linear operators 7~* : A' — > X is called a semigroup on X, if 

(a) — I {I denoting the identity on A'), 

(b) r*+" = r*r" for aii t, s > o. 

Further, if ||T*|| < 1, the family is called a semigroup of contractions. 

If limf_>.o ||7~*/ — /II = for every f E X, then T* is a continuous semigroup (Cq semigroup). 

The transfer operator T'* is a Co semigroup of contractions on L^, see (20j for a proof (in 
particular Remark 7.6.2 for the continuity). 

^/i is absolutely continuous (with respect to m, if not specified further), if there is an < / G L^m) such that 
for all m-measurablc B C M we have fJ,{B) = fg f dm. The function / is called the density of fi. 



V\f {x) = j{'^-'x)\D<^~'xl 




(1) 
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Definition 2.2. For a semigroup T* we define the operator A : 'D{A) — > X by 



with I? (-4) C X being the Unear subspace of X where the above Umit exists. The operator A is 
called the infinitesimal generator of the semigroup. 

For V^, the infinitesimal generator turns out to be (provided the Fi are continuously differen- 
tiable) 

ApFf^-dWifF), (2) 

see [2U]. The following result (see eg. Theorem 2.2.4 [23]) shows the connection between the 
eigenvalues of the semigroup operators and their infinitesimal generator: 

Theorem 2.3 (Spectral mapping theorem). Let T* be a Cq semigroup and let A be its infinitesimal 
generator. Then 

where (t(-) denotes the point spectrum of the operator. The corresponding eigenvectors are identical. 

This has important consequences for invariant densities: 
Corollary 2.4. The function f is a invariant density ofV^ for allt>0 if and only if Appf — 0- 

According to the discussion at the end of Section [2T] we have 
Corollary 2.5. The eigenvalues of App Us on the imaginary axis. 

3 The infinitesimal generator, almost-invariant sets, and es- 
cape rates 

In this section we discuss almost-invariant sets and escape rates via a spectral analysis of the 



infinitesimal generator. As remarked near the end of Section 2.1 the spectrum of V* lies on 
the unit circle, and lacks a spectral gap. Applying the spectral mapping theorem, we see that the 
spectrum of A must be pure imaginary. In the following discussion, to prove formal results, we 
will add a small random perturbation to V^. Later, in the numerics section, we will see that our 
numerical methods introduce a numerical diffusion that plays the role of creating a spectral gap. 

3.1 Stochastic perturbations 

In many real world situations, a deterministic model of some physical system is not appropriate. 
Rather, one should account for the fact that certain external perturbations are present which might 
be unknown or for which a detailed model would be overly complicated. Often, it is appropriate 
to account for these influences by incorporating a small random perturbation into the description. 
From a theoretical point of view, this even facilitates the analysis of the system: Under certain 
assumptions, the transfer operator becomes compact. 

We therefore now leave the deterministic setting and assume that our dynamical system de- 
scribed by the ordinary differential equation dx/dt = F{x) is slightly stochastically perturbed by 
white noise. An exact mathematical treatment of this topic would require tools which are beyond 
the scope of this work (for an exact derivation see [20], chapter 11). We are primarily interested in 
the time evolution of probability densities. The following material highlights the relevant formal 
statements and attempts to point out the intuition behind them. Instead of an ordinary differential 
equation, we now deal with a stochastic differential equation (SDE) 

dX dW 
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with £ > and W being a d-dimensional Wiener process. The solutions are time-dependent 
random variables X(t) with values in M.'^. Just as in the deterministic case, we look at the 
evolution of density functions; now the density functions represent the distribution of random 
outcomes of the variables X{t). The density functions / satisfy 

Prob(X(t) e A) = / f{x,t) dx. 

J A 

Assuming the existence of such a density f{x,t), where /(•, 0) — fo is given, we may again define 
the transfer operator as V^fo :— f{t, •). If the vector field F is smooth enough we have following 
characterization of the density function, the Fokker-Planck or Kolmogorov forward equation, cf. 
[10], 11.6: 

^ = yA/-div(/i^)=: A/. (4) 

Proposition 3.1 ([TJ The operator (with Neumann boundary condition^ is the in- 

finitesimal generator of Cq semigroups Vl i on , Vl 2 on , and Vl q on . 

We note that the operator i is identical to the transfer operator . 

Invariant densities Again, we are particularly interested in distributions (densities) which do 
not change under the evolution. Those are again given by fixed points of Vl i or alternatively, by 
functions in the null space of A^. One can show [3^ that ^ is compact and thus the null space 
of Ae is finite dimensional. Furthermore, due to the white noise, the support of the stochastic 
transition function associated to V^. i is unbounded, and the null space of Ae is one-dimensional, 
i.e. there is a unique invariant density (see again [30j , Theorem 1), characterizing the long term 
dynamical behavior of ([3|. 



3.2 Almost-invariant sets 

We call a set A G M almost-invariant with respect to a (not necessarily invariant) probability 
measure v (cf. [mE]), if 

Ai^y- '^'ZIT^' "' (5) 

for modest times t. The analogous expression for (f> perturbed by a small random perturbation is 

Prob.(X(0) G A) " ^ 

for modest times t. 

We can alternatively characterize this property of a set using an infinitesimal representation 
of Vl- To this end let f £ and consider a measurable set A C M. We define the functional 

Ae,A ■■ ^VAiAe) ^Rhy 

A,Af := lim / ~ ^ dm, f G Va{A,), (7) 

where DAiA^) is the linear subspace of where the above limit exists. Let :— fxA/{J fXAdm). 
Proposition 3.2. Let < / G 2-'^(y^j) be the density of the probability measure v. Then 

pIJA) ^ 1 + {AeMA)t + o{t) (8) 

for measurable A as i —>■ Oj^ 

^The dynamical reason why Neumann boundary conditions are chosen here is discussed later. 
^For two functions /, g : K — >• R we say "/(t) = o{g{t}) as i — ^ 0" , if limt^o = 0. 
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Proof. Wc have 



pIM) - 1 ^ ProW(X(0) e A,Xit) eA)- i^jA) 
t tProb^(X(0) e v4) 

_ L T'IUxa) dm-J^fdm _ f VI f A - Ja 



. .X I dm, 

ti^iA) J A t 

where the second equation follows from the fact that ^(0) is distributed according to v (with 
density /), and ^(0) e A] which gives the density fxA for -'^(O), see Section 3.1 Thus for t — >■ 
the claim follows. □ 

Correspondingly, a set A will be almost invariant with respect to the measure with density 
f ^L\ifA,,AfA-Q- 

There is a strong connection between almost-invariant sets and the spectrum of the generator. 
The following theorem illustrates this with a simple heuristic for producing almost-invariant sets. 

Theorem 3.3. Suppose that the generator possesses a real eigenvalue A < with corresponding 
eigenfunction f . Then J f dm ^ 0. Define A+ = {/ > 0}, A- = {f < 0}. LetfA+,.fA- e T)a{A). 
Then 

Proof. By Proposition 5.7 in [8j and the spectral mapping theorem we have that 

piJA+) + piJA-)^expit\) + l, 



where v is the probability measure with density |/|. Using Proposition 3.2 we obtain 1 + 
Ae,A+ \fA+ \t + I/a- \t = exp(a) + o(t), i.e. 

A,A^ |/^+ I + A,A- \fA-\ = '"'"'^^^^^ ^ + 0(1) 

and for t — > we obtain the claim. □ 
For A ~ Theorem 



3.3 



yields Ai;^a± l/yi± I ~ 0, which means that the sets and A" will be 
almost-invariant with respect to the probability measure with density |/|. Other techniques for 
extracting sets , A~ from the eigenfunction / may be found in [T3 Ej . The papers [13 [T3] 
also discuss almost-invariance with respect to physical invariant probability measures, a property 
that is particularly meaningful when studying typical dynamical behavior. In all cases, the basis 
for these methods are eigenfunctions of A^ corresponding to (real) eigenvalues close to 0. 

If Ae has a complex eigenvalue with real part close to zero, then the corresponding complex 
eigenfunction may also be used to construct almost-invariant sets. 

Lemma 3.4. Let A^f ~ XAf with Xa G C and let fre and /,;,„ denote the real and imaginary 
part of f, respectively. Let t > Q be such that e*^-^ = Xpp G M. Then Vlfre = Xppfre o^nd 

T-'lfim = XpFfim- 



Proof From the proof of Theorem 2.2.4 [23] we have V^f Xppf. Note, that "P* : L^{M,R) O. 

'* f 



By linearity we have V^f = V^fre + i Vlf,^. Thus, 



Xpp fre +i Xppfim — 'Plfre +* "P'/i 



The claim follows immediately. □ 

Hence, if for a i > we have 1 ss e^^^ G M, then the real (and imaginary) part of the 
corresponding eigenfunction yields a decomposition of the phase space into almost invariant sets 
in the sense of Theorem 13.31 
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3.3 Escape rates 

We have seen in the previous section that there is a connection between eigenvalues of Ae close to 
zero and almost-invariant sets. There is also a strong connection between the escape rate of the 
sets A"*", and the corresponding eigenvalue of Ae- 

In the deterministic setting (e — 0) the upper (resp. lower) escape rate of a set A under the 
time-i map of the flow $* is defined as: 

E{A) = -liminf ilogm(An$-*An---n$-'=*A) (9) 
E{A) = -limsup|logm(An$~*74n---n$"''*A) (10) 

fe— >oo k 

If both limits are equal, we call E{A) = \imk_^oo I logm(A n ^~*A n • • • n $~'^M) the escape 
rate from A. The escape rate is the asymptotic exponential rate of loss of m-mass from the set 
A. One might expect that almost-invariant sets have low escape rate and vice-versa, but simple 
counterexamples in HTJ show that this is not the case. This is because the notion of almost- 
invariance as defined above is a finite-time property, while escape rate is an asymptotic quantity. 
The £ > version of Q is 

E,{A) -limmil\ogPmh,n{X{0)e A, X{t)e A,..., X{kt)e A) (11) 

fc— >-oo k 

= -hminf^log / (7'1^)'^(1) dm. (12) 

fc— voo K J 

where ^(/) ■— VK/xa) is the standard restriction of the operator to A [21]. The following 
theorem provides a recipe for constructing sets with low escape rate from eigenfunctions of Ae 
with real eigenvalues. 

Theorem 3.5. Suppose that Aef — A/ for some A < and f £ L°°{m). Then 

Ee{A+) < -tX and E^iA-) < -tX, 
where A+ ^ {f > 0} and A" = {/ < 0}. 

Proof. Let / be scaled such that / |/| dm — 2, and define the (signed) measure i/ as I'iA) — 
f dm. Then = and = are both probability measures. For an event E define 

Probiy(£) := Proby+(£) — Proby-(£). Since / is an eigenfunction of V\ with eigenvalue e'^*, we 
have 



gAfct _ gAfct 



/ fdm = / V^^f dm ^Proh^{X (kt) £ A+) 
Ja+ Ja+ 

= Proh^{X{et)eA+,i^O,...,k) 



fe-i 

^PYoh^{X{nt) e A-,X{et) eA+,i = n+l,...,k) 

n=0^ ^ ' 



It follows for the summands p„ in the above sum: 

p„ = Proh-p^^t^{X{0) e A-,X{et) e A+,£=l,...,k-n) 



e 



Pmh^{X{{)) e A-,x{lt) e ^+,^= l,...,fc-n) 
PvoK-{X{it) e = l,...,A:-n), 



where in the first equation the action of the transfer operator P* on the measure v is defined 
through its action on the density / of v. The second equation follows from v being an eigenmeasure. 
Clearly, p„ is non-positive, and thus 

Prob^(X(^t)e A+,£ = 0,...,fc) >e^'=* 

holds. The rest of the proof follows the lines of the one of Theorem 2.4 in [17]. □ 
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The recipe of Theorem |3.5| is in fact the same as in Theorem |3.3| The difference is the measure 
used: in Theorem|3.3[ almost-invariance is computed with respect to the particular measure v with 



density |/| where / is the eigenfunction in question, whereas in Theorem 3.5 escape is computed 
with respect to Lebesgue measure. 

Remark 3.6. A more natural notion of escape rates in the time-continuous case would be the 
following: 

E^[A) := -liminf-Prob„(X(s) G A, s e [0,t)). 

However, for this definition it is more complicated to obtain similar results to the one in Theo- 
rem 3.5 and a fuller discussion will appear elsewhere. One can view (11 1 with t = 1 as the "time 



sampled version" of this continuous-time definition. 



4 Numerical approximation 

Having seen that certain eigenpairs of the transfer operator (resp. infinitesimal generator) carry 
the information we are seeking, we describe in the sequel our proposed approximation of these 
eigenpairs. To this end, we define finite dimensional approximation spaces and consider the eigen- 
value problem projected onto these spaces. Throughout this section we assume that the underlying 
deterministic vector field F is smooth. 



4.1 Ulam's method 

We describe here the "standard" Ulam approach; see the surveys [13l |6] for more details. We 
partition AI into d-dimensional connected, positive volume subsets {iJi, • . . ,i3ri}- Typically, each 
Bi will be a hyperrectangle or simplex to simplify computations. As an approximation space we 
consider the space A„ = splxsj, . . . ,Xb„} of functions which are piecewise constant on the cells 
of the partition. Let iTn ■ ~> A„, 7r„/ — X]r=i m{B ) Ib f '^"^ XBi, be the L^-orthogonal 
projection onto A„. We let V!^ : A„ — >• A„, '■— 7r„7'*, be the approximate Frobenius-Perron 
operator. Note that V^XBi = T^nV^XBi = ^(gj) Ib, 'P^XBi dm XBj , i-e. the matrix repre- 

sentation G K"^" of with respect to the basis xbi , ■ • ■ , XSn and multiplication on the left 



IS 



iPnh = -7^ / -P'XB^ dm 
JBi 



1 f ,^_ m(i3.n$-*i3,) 
^j/ JB, m{Bj) 



This matrix is easily constructed numerically using eg. GAIO [B]. 

In the stochastic setting, letting „ : A„ — A„, P| „ := T^nVl, one obtains 

as above. In principle one can use Monte-Carlo integration to compute these entries, however, 
this is not particularly efficient. 



4.2 Ulam's method for the generator 

We partition M as in the standard Ulam's method. We will see that the numerical scheme itself 
introduces some diffusion and we therefore consider the deterministic generator A (i.e. e — Q). 

We wish to construct an operator An ■ A„ — >■ A„ that is close in some sense to the operator 
A. Motivated by Ulam's method, one would like to form TTnA, which unfortunately does not exist, 
because A„ ^ 'C'iA), cf. Instead of differentiating w.r.t. time and then doing the projection, 
we swap the order of these operations. Let us build the Ulam approximation first, which will 
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not be a semigroup any more, but for fixed t it approximates P*. Taking the time derivative, our 
candidate approximate operator is 

A^f := lim ( ""^*""/-""^ ) . (13) 

The following lemma emphasizes the intuition behind this definition: if is a suffleiently good 
approximation (for small t) of the Markov jump process generated by the dynamics on the sets 
Bi, then An will be the generator of this process. To be exact, the following is the case: 

Lemma 4.1. The matrix representation of An ■ A„ O with respect to the basis xi, ■ • ■ , Xn under 
multiplication on the left is 

(A ) } t-m{Bj) . .. 

^ "^^^ " ^ miB, n <f>-*i?,) - TO(i3,) , ^^^> 

lim ; — , otherwise. 

t->-0 t ■ Tn{Bi) 

Proof. We consider the action of T'* on ■ 



lim TTn = iim > — , „ , / am y r . 

t^o " t t^o^ m{B.) \Jb. t ' 



lim — ( / dm | y r . 

+ lim — -— - / — - dm — / — - dm y r . 

m{B,n<^>-\Bj)) m(B, n $-*B,) - m(B,) 

> hm — yh. + hm ; — ^ yr. 

^t^o t-m{Bj) t^o t-m{B,) 



Thus under left multiplication we obtain (14). □ 



Remark 4.2. Lemma 4.1 states that {An)ij is the outflow rate of mass from Bi into Bj. 

Let TZl^ :— exp (tAn) — I — T^n + exp(t^„ Ivn)""?! denote the semigroup generated by An. Then 
TZ*^ and are near to each other in the following sense, cf. [T5] : 

Proposition 4.3. As t ^ 

Kf -^nV'f = 0{t^) (15) 

for aZ//e A„0 

The following lemma allows us to construct An without the computation of the flow $* . 

Lemma 4.4. For i ^ j, define to be the the unit normal vector pointing out of Bi into Bj if 
Bi r\ Bj is a d — l-dimensional face, and the zero vector otherwise. The matrix representation of 
An ■ A„ O with respect to the basis xi , . . . , Xn under multiplication on the left is 

— ^— - / ma,x{F{x) ■ nij,0} dmd-i{x), i^j; 
/ . ^ ; m{Bj) Js^nB, 

{An).,j = ( m{B,), ^ , , (16) 

-> — 7T7\\An)ij, otherwise. 



'For two functions /, g : M — >• R we say "f{t) = 0{g(t)) as t — s> 0", if limsupj^Q < oo. 
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Proof. From (|14|) we have for i ^ j that An^ij = Hmt_j.o "^'^^t-iMB -^'^ ■ Denoting Mij{t) = m{Bi n 



<f> we have that An^ij = M[j{Q)/m{Bj) where the prime denotes differentiation with respect 

to t. The quantity M[j{Q) is simply the rate of flux out of Bi through the face Bi n Bj into Bj 
and so M^(0) = /g.^^^ max{F(a;) • nij,Q)dmd-i{x). 

For the diagonal elements A„ we have AnM — limt_>.o 



m(Bin^~* Bi)-m(Bi) 
t-m(Bi 



Note that m{Bi 



m{B, n <^-'Bi) - m{B, \ <S>-'B,). Clearly B, \ = B, D U,_,, = U,^.^. n 



modulo sets of Lebesgue measure zero. Thus, m{Bi) — ra{Bi n $ 

5^m(S,n<I.-'S,) 



Now, by (14 1, A„,ii = -lim* 



m{Bi) 



J^i m(Bi)^n,ij- 



An,' 



In one dimension, ( 16 1 has a particularly simple form. 



Corollary 4.5. Let M = . Assume (without los^ that F{x) > 0. Let = xq < xi < 

Xn — 1 and consider the partition {-Bi, . . . , -B„} with Bi — [xi„i, Xi\. Then 



□ 



< 



F{x^)/miBj), j = i + l; 
{An)ij = { -F{xi)/m{Bi), j = i; 

0, otherwise 



(17) 



Remark 4.6 (Connections with the upwind scheme). Clearly, An is the spatial discretization 
from the so-called upwind scheme in finite volume methods; cf. pij . The scheme is known to be 
stable. Stability of finite volume schemes is often related to "numerical diffusion" in them. 

We demonstrate numerical diffusion on a simple example, for further details we refer to |21j 
Section 8.6.1. Set M = and F{x) = F > Va; e M. Then Af = -Ff and An is the backward 
difference scheme. Consider /„ :~ TTnf as a vector of values. For sufficiently smooth / we have 



{Anf)^ = nP - fn,^) = (vT^^/)^ + O (^-1) 



but 



{Anf)i = nP [fn.i-l - fn.i) 



nF 
P 



fn,i—l fn.i-\-l , fn,i—l ^fn.i ~t~ J] 




n.i-\-l 



2 2 

fn.i—1 fn.i-\-l , ^ fn,i — l ^ ^fn.i ~^ /"n/i+l 

2nr^i 



hence Anf is a better approximation of A^f, with % = and Ae is the infinitesimal generator 
associated with the SDE ^ = Fi^) + £^jr than of Af. That is why one expects quantities 
computed by An to reflect the actual behavior of A^ . Since the above equations contain only local 

estimates, for a general nonconstant F is the diffusion also spatially varying; i.e. '^^^ = ^i^- 
Note, that this notion of numerical diffusion is easily extended to multiple space dimensions. 



The definition ( 13 ) allows us to interpret the numerical diffusion in yet another way. We 
showed in Proposition 4.3 that is the transition matrix of a Markov process near the Markov 



jump process generated by An for small t > 0. The discretized FPO P,* can be related to a non- 
deterministic dynamical system, which, after mapping the initial point, adds some uncertainty to 
produce a uniform distribution of the image point in the box where it landed; see |12|. Thus, this 
uncertainty resulting from the numerical discretization, equivalent to the numerical diffusion in 
the upwind scheme, can be viewed as the reason for robust behavior — stability. 



Finally, we show that our constructions (16) and (17) always provide a solution to the system 
Anf = for some / G A„. 



^If _F ^ and -F ^ 0, we have one or more stable fixed points, and every trajectory converges to one of them. 
Hence, there is no interesting statistical behavior to analyze. 
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Lemma 4.7. There exists a nonnegative, nonzero f G A„ so that Anf = 0. 
Proof. Let Dn^ij = m{Bi)5ij and note that '■= D^^ A^Dn satisfies 



{Qnh^l "''^'"ZiB., , , . . (18) 



m{B,) * 7^ 

- Ej5^» Sl7y(^")«J' Otherwise. 



Note that all row sums of Qn equal zero. Let c ~ '^j^iQn,ij- The matrix Q„ := Qn + c/ 
is nonnegative with all row sums equal to c. By the Perron-Frobenius theorem [3], the largest 
eigenvalue of Qn is c (of multiplicit} ^ possibly greater than 1) and one of the corresponding left 
eigenvectors pair u„ is nonnegative. Clearly m„ is a left eigenvector of Qn corresponding to the 
eigenvalue and thus DnUn is a nonnegative left eigenvector corresponding to for An- □ 

4.3 Spectral collocation for the generator 

The eigenfunctions of A^ are smooth. This motivates the use of smooth approximation functions, 
e.g. polynomials. We here outline the general principles of spectral collocation methods; for a 
more thorough presentation we refer to [5], [3] or |27) . 

Choose a family of approximation spaces {V^}„gN, such that Vn C C°°{M) for all n. Depending 
on the type of the phase space, we use two different approximation spaces. We introduce them in 
one dimension; the multidimensional ones can then be constructed by tensor products. In both 
cases, the approximation space comes with an associated set of collocation nodes: 

• Periodic domain/uniform grid. We have M = and restrict ourselves to odd values of n. 
Then the basis we choose for Vn is 

\ / -n/2-l<fc<n/2 ■ 

The associated collocation nodes are the uniform grid {0, 1/n, . . . , (n — l)/n}. 

• Standard domain/Chebyshev grid. Here, AI = [—1,1]. The space Vn is spanned by the 
monomials of order to n. We use Chebyshev polynomials as basis functions: 

{cos (k arccos (a;))}o<^<„ , 

together with the Chebyshev grid {— cos{2TTj /n), j = 0, . . . ,n}, as collocation nodes. 

Let f £ Vn and Z„ : C°° — Vn be the interpolation operator for the given collocation nodes. We 
define the approximate generator by 

A^_nf • — -^nA^ f . 

For both cases we have following: 

Theorem 4.8 (Spectral accuracy, [5]). For f £ C°°{M) let fn be the best approximation of f in 
Vn w.r.t. the supremum norm \\ ■ ||oo- Then for each fc G N there is a Ck > such that 

||/-/„||oo <cfc forallneN. (19) 

Convergence then follows from standard results on the analysis of Galerkin methods for elliptic 
differential operators and spectral approximation (cf. also jl8j): 

Theorem 4.9. Let M be a compact tensor product domain with smooth boundary. Let F be a 
smooth vector field on M andAe, Ae.n defined as above, with Neumann (resp. periodic) boundary 
conditions. Then the eigenvalues and eigenfunctions of Ae are approximated by the corresponding 
eigenvalues and eigenfunctions of Ae^n with spectral accuracy. 



®If Qn is aperiodic (there exists k such that > 0) then the eigenvalue c is simple and the corresponding 
eigenvector is positive (see [3]). 
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Why Neumann boundary conditions? The objects that we are approximating are densities, 
so there should be no "loss of mass" as time goes on. In a closed physical system, the flow normal 
to the boundary is zero. Hence, there is no loss of mass through advection. The physical meaning 
of diffusion is a flow of mass in the direction opposite to the gradient and proportional to its 
magnitude. Thus, no loss of mass via diffusion translates into Neumann boundary conditions for 
the densities: their normal gradients have to be zero at the boundary. Equivalently, one may write 
Q as a continuity equation, 

dj = div(7), 

with I :— ^V/ — fF being the probability flow or probability current. The condition I\qm = 
leads precisely to Neumann boundary conditions. 



4.4 Algorithms 



Algorithm 1 (Ulam's method). 

1. Partition M into connected sets {Bi, . . . , i?„} of positive volume. Typically, each Bi will be 
a hyperrectangle or simplex. 



2. Choose t > and compute the matrix P"* or P|„ as described in Section 



4.1 



3. Estimates of invariant densities for $* are given by fixed points of (or P| „): Let uP^ = v. 
Then / ViXB, satisfies V^f = /. 

4. Similarly, eigenvectors of P* corresponding to real eigenvalues A « 1 provide information 
about almost-invariant sets (cf. Theorem 3.3 1 and sets with low escape rates (cf. Theorem 



3.5[ ). 

Note that P^ typically is a sparse matrix and the number of nonzero entries per column will 
be determined by Lipschitz constants of 

Algorithm 2 (Ulam's method for the generator). 

1. Partition M into connected sets {Pi, . . . , P„} of positive volume. Typically each Bi will be 
a hyperrectangle or simplex. 

2. Compute 

— — / max{P(x) • nij,0} dTO<j_i(a:;), i ^ j, 
l^jj JB,nB, 



MB]) JB^nB, 
m{Bi) 



EmlBj) , ^ , , . 

(R \ otherwise. 



where one uses standard quadrature rules to estimate the integral. 

3. Estimates of invariant densities lie in the null space of An- Let vAn — (the existence of 
such a w is guaranteed by Lemma 4.71, then / := X^ILi ''^iXBi satisfies Anf — 0. 

4. Similarly, eigenvectors of An corresponding to large real eigenvalues A < provide informa- 
tion about almost-invariant sets and sets with low escape rates. 

Note that the discretized generator A„ is a sparse matrix since An.ij = if Bi and Bj do not 
share a boundary. 

Algorithm 3 (Spectral collocation for the generator). 
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1. Set up a grid xq, ■ ■ ■ , a;„, which is the Chebyshev grid Xj = — cos (27rj/n), j = 0, . . . ,n, if 
M = [— 1, 1], or the equispaced grid Xj = j/n, if M = T^. (If the state space M is an affine 
transformation of the above ones, the grid is transformed analogously.) For multidimensional 
tensor product spaces the grid is the tensor product of the one dimensional grids. 

2. Let i£ Q , . . . , in} denote the Lagrange basis on the grid, i.e. all £i are polynomials of maximal 
degree n with £i{xj) — Sij. Denote the interpolation on the grid by I„, and define the 
discretization matrix obtained by collocation as 

Note, that the computation of these matrix entries is very simple. If M = T-'^, we may switch 
between the evaluation space and the frequency space simply by the fast Fourier transfor- 
mation. Multiplication by the vector field F is pointwise multiplication in the evaluation 
space, computing the derivative is a diagonal scaling in the frequency space. This can be 
simply extended to M = [0, 1] (cf. [27] Chapter 8), as well as for a multidimensional phase 
space. 

3. Compute the left eigenvector v of ^£,„ at the eigenvalue of smallest magnitude. Then, 
S"=o approximates the invariant density. 

4. Similarly, eigenvectors of A^.n corresponding to large real eigenvalues A < provide infor- 
mation about almost-invariant sets. 



4.5 Solving the eigenproblem 

Once the matrix approximation of the operator or the generator has been computed, we then have 
to solve the corresponding eigenvalue problem. In Ulam's method one seeks dominant eigenvalues, 
and standard Arnoldi type methods easily provide the interesting part of the spectrum (we use 
the eigs function in Matlab). 

For the approximate generator, we look for eigenvalues with small modulus. Here, the Arnoldi 
iteration in eigs uses inverse iteration, i.e. repeatedly solves linear systems of the type A„u = b. 
Since the matrix An is singular, this is not a well posed problem and we might expect some diffi- 
culties with these iterative methods. A possibility for overcoming them is to solve the eigenvalue 
problem for An '■= An — al„xn, with a suitable shift a. This merely introduces a shift of the 
spectrum and the eigenfunctions stay the same. For n big enouglj^ the spectrum of An is going 
to be well separated from zero. The size n of the linear system may be large, so a direct solution 
via LU factorization may not be feasible. However, in Ulam's method for the generator. An is a 
sparse matrix, hence iterative methods (e.g. GMRES) can be used. 

For spectral collocation, in general, the matrix A„ will not be sparse. Nevertheless, in the 
case where the domain is periodic, the matrix-vector product AnX can be computed fast using the 
FFT. 



4.6 Computational complexity 



For d dimensional systems in general, ( |16[ ) shows us that setting up the approximate generator 
requires the computation of (d — 1)- dimensional integrals, while to set up the Ulam-matrix, we 
compute d-dimensional integrals. Of course, the gain is biggest in one dimension (zero dimensional 



integrals are one-point-evaluations). Equally importantly, (16) does not require integration of the 
vector field. 

For Ulam's method and Ulam's method for the generator, the matrix entries should be com- 
puted to an accuracy 0{n^^^'^); otherwise we cannot expect the approximate eigenfunctions to 



An is obtained from the Ulam-based approach, the spectrum of An, which is the generator of a Markov jump 
process, is in the left complex half plane. This does not hold in general for the spectral collocation approach, here 
we rely on the fast convergence of eigenvalues. 
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exploit the full potential of the approximation space (which is an error of 0{n~^/'^) for n partition 
elements in a d dimensional space). Assuming the quadrature rule used to compute the entries to 
also suffer from the curse of dimension, we expect to need 0{n) and ©(n^''^^'/'^') function evalu- 
ations for each matrix entry in the Ulam-matrix, and the matrix arising from Ulam's method for 
the generator, respectively. This explains the third column of Table [l] below. 

5 Examples 

A collection of examples is presented in the following. The examples range over three phase 
space dimensions, different boundary types, and chaotic and regular dynamics. We find that 
the generator approaches are typically superior to the standard Ulam approach, producing more 
accurate results with less computation time. If the vector field is infinitely differentiable, the 
spectral collocation approach can be extremely accurate. In our final example, we find that for 
fiows on complicated attractors, the standard Ulam method can take advantage of the attracting 
dynamics to produce better estimates of invariant densities than the generator approaches. In all 
cases, the operator based methods are far superior to direct trajectory integration for estimating 
the invariant density. Other structures such as almost-invariant sets can only be identified using 
the operator approaches. 

5.1 A flow on the circle 

We start with a one dimensional example, a flow on the unit circle. The vector field is given by 

F{x) = sin(47rx) + 1.1, 

X G = [0, 1] with periodic boundary conditions, and we wish to compute the invariant density 
of the system. Recall that an invariant density / G L^{T^) needs to satisfy Af — 0, where 
Af{x) — — d{fF){x)/ dx. The miique solution to this equation is f*{x) — C/F{x), C being a 
normalizing constant (i.e. such that ||/||i — 1). The honours thesis "SB numerically investigated 
the estimation of the invariant density for this flow, and other two-dimensional flows, by finding 
approximate eigensolutions of the infinitesimal generator using finite difference and finite element 
methods. Here we use the three methods discussed in the previous section to approximate /* and 
compare their efficiency relative to one another and a histogram of a long simulation: 

1. the classical method of Ulam for the Frobenius-Perron operator, 

2. Ulam's method for the generator and 

3. spectral collocation for the generator. 

As we have a periodic domain and the vector field is infinitely smooth, we expect spectral colloca- 
tion to perform very well. Figure [l] shows the true invariant density (dashed line), together with 
its approximations by the four methods. 

5.1.1 Matlab code 

In this ID case, the three methods can each be realized in a few lines of Matlab. For illustration 
purposes, we include the code here. 

Ulam's method 

F = 8(x) sin(4*pi*x) + l . 1 ; definition of the vector field 

n = 32; x = l/(2*n"2) : l/(n"2) : l-l/(2*n"2) ; '/, nodes for spatial integration 

[t,y] = ode23(a(t,x) F(x) , [0 2/n],x'); '/. time integration of the nodes 

I = ceil (max(n*mod(y (end, :), 1) , 1) ) ; '/, cell indices of image points 

J = reshape (ones (n, 1) * (1 :n) , 1 ,n*n) ; construction of transition matrix 
P = sparse (I , J , l/n,n,n) ; 
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[v,d] = eig(f ull (P) ) ; '/, spectrum of transition matrix 

f_n = @(x) abs (v(ceil (inax(x*n, 1) ) , 1) ) ; plot approximate inv. density 

LI = quadl(f_n,0,l) ; f plot (0 (x) f _n(x) /LI , [0 , 1] ) 

The error in Ulam's method decreases Uke 0{n^^) for smooth invariant densities [10 . Thus, 
we need to compute the transition rates between the intervals to an accuracy of 0{n^^) (since 
otherwise we cannot expect the approximate density to have a smaller error) . To this end, we use 
a uniform grid of n sample points in each interval. This leads to 0{n^) evaluations of the vector 
field. For the numbers in Figure [T] we only counted each point once, i.e. we neglected the fact that 
for the time integration we have to perform several time steps per point. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 1: True invariant density (dashed line), approximation by histogramming a long simulation 
(top left, 1024 iterates), Ulam's method (top right, 32 sample points per partition element, 32 
partition elements), Ulam's method for the generator (bottom left, 32 partition elements) and 
spectral collocation (bottom right, 32 collocation points). 



Ulam's method for the generator 

F = aCx) sin(4*pi*x)+l.l; 
n = 32; Fx = F(0: 1/n: 1-i/n) ; 

A = n*sparse(l :n, 1 :n,-Fx) + sparsed :n, [2:n 1] ,Fx) ; 
[v,d] = eigCfulKA')) ; 
[d,I] = sort(diag(d)) ; 

f_n = (a(x) abs(v(ceil(max(x*n,l)) ,1(1))) ; 

LI = quadl(f_n,0,l) ; f plot (0 (x)f _n(x) /LI , [0 , 1] ) 



definition of the vector field 
"it, evaluation on the boundary nodes 

assembling the discrete generator 
"it, spectrum 

"it, plot approximate inv. density 



Here, only one evaluation of the vector field per interval is needed as in one dimension the 
integration on the boundaries reduces to the evaluation of a single boundary point. On a partition 
with n intervals, this method then yields an accuracy of 0(n~^). Note that from Corollary 4.5 it 
follows that the vector v with Vi — 1/F{xi) is a left eigenvector of the transition matrix (17) for 
the generator at the eigenvalue 0. This fact proves pointwise convergence of the invariant density 
of the discretization towards the real one. Thus Ulam's method for the generator is very accurate 
and stable; the error is 0{l/n) as this is the rate at which approximants created from a basis 
of characteristic functions converge in to regular functions. 
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Spectral collocation 



F = a(x) sin(4*pi*x)+l.l; 




% definition of the vector field 


n = 32; Fx = F(0: 1/n: 1-1/n) ; 




'/. evaluation in collocation nodes 


E = ifft(eye(n)) ; 




'/, basis 


FE = fft(diag(Fx)*E) ; 




'/o multiplication by vector field 


I = [0:ii/2-l -n/2:-l] ; 




'/o frequencies 


D = (2i*pl)*I' ; 




% differentiation matrix 


A = -D*ones(l,n) .*FE; 




% discrete generator in frequency space 


[v, lambda] = eig(A); 




'/o spectrum 


f_n = (a(x) real(exp(x'*I*2i*pi)*v(: 


, end) ) ; 


"h approximate density 


LI = quadl(f_n,0,l) ; f plot (a (x) f _n( 


x) /LI, [0,1]) 


■/, plot 



Here, the vector field is evaluated once per grid point (cf. the second line of the code), 
predicted by Theorem 4.8 the accuracy increases exponentially with n (cf. Figure l2|. 



As 



5.1.2 Computational efficiency 

In Figure [2] (left) we compare the efficiency of the four methods in terms of how the L^-crror of 
the computed invariant density depends on the number of evaluations of the vector field. The 
L^-error was computed using an adaptive Lobatto quadrature as implemented in, e.g., Matlab's 
quadl command. These errors are due to the approximation space, the approximate numerical 
computation of matrix entries, and due to solving the discretized eigenvalue problem. To illus- 
trate how these effects accumulate, in Figure [2] (right) we compare the errors of three different 
approximations, for different (uniform) partitions: the approximate invariant densities obtained 
from the two Ulam type methods, and the projection of the true invariant density. We conclude 
that the error due to the approximation space dominates the total error. 




# of evaluations # of partition elements 



Figure 2: Left: L^-error of the approximate invariant density as a function of the number of 
evaluations of the vector field (2'* - 2^ partition elements, resp. collocation points, 2® - 2^^ number 
of iterations). Right: L^-error of the approximate invariant densities obtained by projection of 
the true invariant density on the underlying partition, by Ulam's method, and by Ulam's method 
for the generator. 

The number of grid sets (32) chosen in Figure [l] is tiny, and chosen merely for illustration 
purposes. Likewise, the maximum number of evaluations used in both generator schemes of 512 
is also obviously tiny, and in practice one could very cheaply increase the number of grid sets and 
concomitantly the number of evaluations. The main message from this example is that in the right 
setting: low dimension, periodic domain, and infinitely smooth vector field, spectral collocation 
can significantly outperform standard Ulam and generator Ulam. Ulam's method for the generator 
is clearly outperforming standard Ulam in this example, and we will see that it continues to do so 
across a variety of dimensions, domains, and vector fields. 
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Table 1: Sources of computational cost. The dimension of the approximation space is denoted 
by n, d denotes the dimension of state space. In contrast to the other two methods, spectral 
collocation produces a full matrix. One can reduce the cost for solving the eigenproblem in this 
case to e'(nlog(n)) • #(GMRES iterations) by using an iterative solver (e.g. GMRES), and FFT 
to compute matrix- vector products. 



method 


approximation 
error 


flops to set up matrix / 
time integration 


dimension 
of spatial 
integrals 


flops to solve EVP / 
vector iteration type 


Ulam's method 




/ yes 


d 


0{n) 1 forw. 


Ulain's method 
for the generator 




0(„(2d-l)/d) / 


d-\ 


0{n) 1 backw. 


spectral coUocation 
for the generator 




0{n) 1 no 





Oijv') 1 backw. 



5.2 An area-preserving cylinder flow 

We consider an area-preserving flow on the cylinder, defined by interpolating a numerically given 
vector field as shown in Figure [sj which is a snapshot from a quasi-geostrophic flow, cf. [2S] . The 
domain is periodic with respect to the x coordinate and the field is zero at the boundaries y = 
and y — 10^. Again, we apply the three methods discussed in Section |4] in order to compute 
approximate eigenf unctions of the transfer operator and the generator. 




Figure 3: Vector field of the area-preserving cylinder flow. Shown are also two saddle equilibria 
(black dots) together with their stable and unstable manifolds. 

5.2.1 Perturbing the model 

Since the flow is area-preserving, we have a continuum of closed orbits, and each of them is an 
invariant set. On each, an invariant measure can be supported, i.e. the eigenspace at 1 (resp. 0) 
is infinite-dimensional. Both Ulam type methods can be interpreted as a small random perturba- 
tion of the original system (albeit different ones), yielding a unique stationary eigenvector of the 
resulting matrix. However, simultaneously with decreasing the cell size, the size of this "built-in" 
perturbation shrinks and it is therefore not immediate to which invariant measure the discretiza- 
tion schemes converge to. In order to make the problem well posed and the results comparable 
among the three methods we therefore need to artificially add a random perturbation (in form of 
a diffusion term) to the model. 

We therefore choose the noise level e such that the resulting diffusion coefficient has 
the same order of magnitude as the numerical diffusion present within Ulam's method for the 
generator, but is still larger than the latter. Since the estimate from Remark [4 . 6| yields a numerical 
diffusion coefficient of « 120 (see also [H]), we choose £ = \/2 ■ 500 here. 

In Ulam's method, for the simulation of the SDE (|3| a fourth order Runge-Kutta method 
is used, where in every time step a properly scaled (by a factor \/h ■ e, where h is the time 
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step) normally distributed random number is added. We use 1000 sample points per box and the 
integration time T = 5 • 10^, which is realized by 20 steps of the Runge-Kutta method. Note that 
the integrator does not know that the flow lines should not cross the lower and upper boundaries 
of the state space. Points that leave phase space are projected back along the y axis into the 
next boundary box. An adaptive step-size control could resolve this problem, however at the 
cost of even more right hand side evaluations. The domain is partitioned into 128 x 128 boxes. 
For Ulam's method for the generator, again, we employ a partition of 128 x 128 boxes and 
approximate the edge integrals by the trapezoidal rule using three nodes (which turns out to be 
sufficiently accurate in this case). For spectral collocation, we employ 51 Fourier modes in 
the X coordinate (periodic boundary conditions) and the first 51 Chebyshev polynomials in the 
y coordinate, together with Neumann boundary conditions. Due to spectral convergence, this 
smaller number of modes already yields an accuracy in the invariant density comparable to the 
one of the other two approaches. 

5.2.2 Computing almost invariant sets 

In Table |2] we compare approximations for the three leading eigenvalues from our three schemes. If 
we assume that spectral collocation gives a very accurate result, Ulam's method should be 0{'nr^) 
from it. With n = 128 this matches well with the numbers in the table. Ulam's method for the 
generator generates additional numerical diffusion; this scales the corresponding eigenvalues. In 



Table 2: Approximate eigenvalues. 



method 


A2 


A3 


A4 


Ulam's method (log(A,)/T) 
Ulam's method for the generator 
spectral collocation for the generator 


-1.64- 10-** 
-1.98- 10-** 
-1.65 • 10-** 


-0.91 • lO--" 
-1.03-10-^ 
-0.92 • 10-^ 


-1.06- 10-^ 
-1.19- 10-^ 
-1.05 • 10-^ 



Figure [4] we compare the approximate eigenvectors of the second, third and fourth relevant eigen- 
values of the transfer operator (resp. generator). Clearly, they all give the same qualitative picture. 
Yet, the number of evaluations of the vector field (and thus the associated cpu times) diff'er sig- 
nificantly, as shown in Table [s] The regions enclosed by the stable/unstable manifolds of the two 

Table 3: Number of evaluations of the vector field and cpu times (on a 2.5 GHz Core Duo) for 
computing the approximate operator resp. generator and cpu times for computing the leading 4 
eigenvalues /-vectors. 



method 


^ of rhs evals 


time matrix 


time eigs 


Ulam's method 


w 3 • 10« 


3020 sec. 


1.0 sec. 


Ulam's method for the generator 


w 3 • 10^ 


143 sec. 


0.8 sec. 


spectral collocation for the generator 


w 3 • 10^ 


0.4 sec. 


3 sec. 



saddle equilibria shown in Figure [3] are obviously invariant for the noiseless flow. The introduction 
of noise means that it is possible (but with a rather low probability) for stochastic trajectories to 
jump between these regions, making these regions almost-invariant sets (cf. Section [3.2[ ). Figure 
|4] shows that these almost-invariant sets are highlighted in the 2nd, 3rd, and 4th eigenvectors as 
level sets at high gradient regions in the eigenvectors (in this example, near the yellow bands) . The 
2nd and 3rd eigenvectors highlight regions near different stable/unstable manifolds at the yellow 
bands and the 4th eigenvector shown in Figure [s] shows a combination. We refer the reader to 
[9lll6j for details on the connection between invariant manifolds and almost invariant sets. As one 
sees in Figure [5j the eigenvectors from Ulam's method for the generator are smoother compared 
to the ones from the standard Ulam method. Ulam's method for the generator appears to be more 
accurate, at least for this system with a small amount of noise. A possible explanation is that 
with Ulam's method for the generator we formally average the diffusion, while with the standard 
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(a) Ulam's method 




(b) Ulam's method for the generator 




(c) Spectral collocation for the generator 



Figure 4: Quasi-geostrophic flow: Eigenvectors at the second, third and fourth eigenvalue (from 
left to right). 
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Figure 5: Comparison of the approximate eigenvector at A4 by Ulam's method (left) with the one 
computed by Ulam's method for the generator (right, together with the two saddle equilibria and 
their invariant manifolds). 

Ulam method we simulate a finite number of random trajectories, which incompletely sample the 
true noise distribution. 

5.3 A volume-preserving 3D example: the ABC-flow 

In this section we demonstrate the efficacy of our generator approach for a three-dimensional flow. 
We consider the so-called ABC-flow p], given by 

X = asin(27rz) + ccos(27ry) 
y — 6sin(27ra;) + acos(27rz) 
z = csin(27ri/) -f 6cos(27ra;), 

on the 3 dimensional torus. The flow is volume-preserving, so in order to make the computation 
of the spectrum well posed we need to add some artifical diffusion again. The vector field is very 
smooth, and in numerical experiments a diffusion coefhcient e ~ 0.04 turned out to work well. 
For a = \/3, h = yjl and c — \, there is numerical evidence that the flow exhibits complicated 
dynamics and invariant sets of complicated geometry [TTl [TB] . 

We approximate the leading few eigenfunctions of Ae with Ulam's method for the generator 
on a 64 X 64 X 64 grid. We used Gauss quadrature with 16 nodes on each face of a box in order 
to approximate the surface integrals, requiring ~ 16 • 64^ « 4 • 10^ evaluations of the vector field. 
The L^-error of the approximate invariant density is « 10~^^. The second and fifth eigenfunctions 
are shown in Figure |6] The regions highlighted in red and blue are invariant cylinders of the 
deterministic flow; see |16j for more details on how to extract invariant and almost-invariant sets 
from the eigenfunctions. 

We also apply spectral collocation onallxllxll grid, requiring ll'^ = 1331 evaluations 
of the vector field. The L^-error of the approximate invariant density is « 10~^^. We observe 
fast convergence (after comparing with computations on coarser resolutions) of not just the first 
nontrivial eigenfunctions, but also of higher order ones. 
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Figure 6: Second and fifth eigenvectors of the infinitesimal generator - Ulam type discretization. 
Left: shces, right: regions, where the absolute value of the left eigenfunctions exceed a given 
threshold - the "core" of the almost invariant sets, representing invariant cylinders in this example. 

5.4 A dissipative 3D example with complicated geometry: the Lorenz 
system 

As a last example we consider a system where the effective dynamics is supported on a set of 
complicated geometry and fractional dimension, the well-known Lorenz system 

i = (T(y - x), 
y = x{g - z) - y, 
z = xy — j3z, 

where we make the standard parameter choices cr = 10, = 28 and /3 = 8/3. Since the effective 
(complicated) dynamics happens on a lower dimensional attractor we do not expect spectral 
collocation to work well and we therefore apply Ulam's method for the generator here. 

A decade ago, numerical techniques were introduced for computing box coverings for attractors 
of complicated structure, cf. [71 [S] . These techniques make use of the fact that the set M to be 
computed is an attractor, hence each trajectory starting in its vicinity will be pulled to M in a 
fairly short time. In our approach time is not considered, we only use the vector field. Since 
the boundary of the box covering does not have to coincide with the boundary of the attractor, 
a tight box covering might not show the desired results, because of relatively big outflow rates 
in boundary boxes. The simplest idea is to use a big enough rectangle - in our case [—30, 30] x 
[—30, 30] X [—10, 70]. However, this set is not invariant, hence the question arises, which boundary 
conditions one should impose. Translating the condition "no probability flow on the boundary" in 
this Ulam type approach leads to the condition that the flow rates at the boundary are ignored. 
The attractor is then extracted by simple thresholding: regions where the approximate invariant 
density is nonzero are expected to lie in a small neighbourhood of the attractor. The finer the 
resolution, the smaller the diffusion introduced by the discretization and thus the tighter this 
neighbourhood of the attractor is. Having constructed the attractor, we may restrict the other 
eigenfunctions to this set in order to obtain almost invariant sets within the attractor. 
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Figure 7: Second and fifth eigenfunction of the generator computed by spectral collocation. Left: 
slices, right: regions, where the absolute value of the left eigenfunctions exceed a given threshold 
- the "core" of the almost invariant sets. 

For our computation, we employed a grid of 128 x 128 x 128 boxes and used Gauss quadrature 
in 5 X 5 nodes on the faces in order to set up the approximate generator. The attractor is extracted 
by thresholding the approximate invariant density u: We cut off at a threshold value 5 • 10~^ such 
that 96% of the invariant measure is supported on {ui > c}. Figure [s] shows the corresponding 
covering as well as the sign structure of the second and third eigenvector. 




Figure 8: Ulam's method for the generator in an example with a low-dimensional attractor: The 
Lorenz attractor (left) and almost invariant sets: sign structure of the second eigenvector (middle) 
and third eigenvector (right). 



6 Conclusion 

We introduced the infinitesimal generator as a tool for studying long term behavior of flows without 
requiring any trajectory integration. The (typically unit dimensional) null space of the generator 
is spanned by the invariant density of the flow; the invariant density describes the asymptotic 
distribution of trajectories. We showed that eigenfunctions of the generator corresponding to 
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eigenvalues close to had strong connections with almost- invariant or metastable behavior, and 
could be used to identify regions in the flow from which the escape rate is very low. 

We proposed two new numerical approaches for the estimation of the generator; one based on 
Ulam's method and the other based on spectral collocation. We tested our numerical approaches 
on "full phase space" flows in one, two, and three dimensions and found that both approaches 
signiflcantly outperformed the standard Ulam approach based on the Perron-Probenius operator. 
Moreover, all operator/generator based methods were superior to standard trajectory integra- 
tion for estimating the invariant density; the identification of almost-invariant sets is essentially 
impossible without an operator/generator approach. The spectral collocation approach worked 
particularly well in periodic domains when the vector field is infinitely smooth. 
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